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Abstract 

We study the effects of superconducting pairing in small metallic grains. We show that in the 
limit of large Thouless conductance one can explicitly determine the low energy spectrum of the 
problem as an expansion in the inverse number of electrons on the grain. The expansion is based 
on the formal exact solution of the Richardson model. We use this expansion to calculate finite 
size corrections to the ground state energy, Matveev-Larkin parameter, and excitation energies. 



I. INTRODUCTION 

Since mid 1990's, when Ralph, Black, and Tinkham succeeded in resolving the discrete 
excitation spectrum of nanoscale superconducting metallic grains^, there has been consid- 
erable effort to describe theoretically superconducting correlations in such grains (see e.g. 
Ref. |2j for a review). A key question in any such description is how results of the BCS 
theory are modified in finite systems. In this paper we address this problem by developing 
a systematic expansion in the inverse number of electrons on the grain for the low energy 
spectrum of the problem. 

In the absence of spin-orbit interaction and magnetic fields one can describe^ supercon- 
ducting correlations in weakly disordered grains by a simple pairing (BCS) Hamiltonian 

n 

Hbcs = Y, e ^W - Xd ^ 4i c h c rt c ji (i-i) 

i,cr i,j = l 

Here e^ are orbital energy levels and d is the mean level spacing, d = (e$+i — e;). Operators 
c ia (cio) create (annihilate) an electron of spin projection a in orbital state i, n is the total 
number of levels, and A denotes a dimensionless coupling constant. The interaction part of 
Hamiltonian (jl.ljl allows only transitions of singlet electron pairs between the orbitals. 

BCS Hamiltonian (jl.lj) is known to be integrable^ and solvable by Bethe's Ansatz. The 
exact solution 5 - yields a complicated set of coupled polynomial equations (see Eq. (J1.2JI 
bellow). As a consequence, very few explicit results have been derived and most studies 
resorted to numerics^^ based on the exact solution. The purpose of the present paper is 
to remedy this situation in the regime when the level spacing is the smallest energy scale in 
the problem. 

BCS Hamiltonian (jl.lj) was studied extensively in 1960 's in the context of pair corre- 
lations in nuclear matter 8 . A straightforward but important observation was that singly 
occupied orbitals do not participate in pair scattering^. Hence, the labels of these orbitals 
are good quantum numbers and their contribution to the total energy is only through the 
kinetic term in BCS Hamiltonian (jl.lj) . Due to this "blocking effect" the problem of di- 
agonalizing Hamiltonian (jl.lj) reduces to the subspace of orbitals that are either empty or 
doubly occupied - "unblocked" orbitals. The latter problem turns out to be solvable 5 by 
Bethe's Ansatz. The spectrum is obtained from the following set of algebraic equations for 
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where m is the total number of singlet pairs and n now is the number of unblocked orbitals €k- 
Bethe's Ansatz equations (jl.2j) for BCS Hamiltonian (jl.lj) are commonly referred to as 
Richardson's equations. Eigenvalues of BCS Hamiltonian (jl.l|) are related to Richardson 
parameters, Ei, via 

m 

E = 2Y,E t + Y<£B (1.3) 

i=l B 

where J2b € b is a sum over singly occupied - "blocked" orbitals. 

In 1977, Richardson used exact solution ()1.2|) to outline^ a method for expanding the low 
energy spectrum in powers of the inverse number of pairs, 1/m. Richardson showed that BCS 
results^ for the energy gap, condensation energy, excitation spectrum etc. are recovered 
from exact solution (jl.2j) in the thermodynamical limit. The proper limit is obtained by 
taking the number of levels, n, to infinity, so that nd — > 2D = const, m = n/2, where D is 
an ultraviolet cutoff usually identified with Debye energy. In particular, for equally spaced 
levels €i, the energy gap A and the ground state energy in the thermodynamical limit are 

Ao(A)= . f ^ s C5 (A) = - J Dmcothl/A (1.4) 

smn(l/A) y 

In the present paper we show that the ground state and excitation energies of BCS 
Hamiltonian (jl.lj) can be evaluated explicitly to any order in d/Ao ~ 1/m in terms of 
the BCS gap Ao, chemical potential /i, mean level spacing d, ultraviolet cutoff D, and the 
thermodynamic density of states v{e). In the physical limit A /D —>■ 0, the expansion is 
applicable for Ao > d. In fact, we believe that in this limit the expansion is in powers of 
d/Ao with a convergence radius dj Ao ~ 1. 

BCS Hamiltonian (jl.lj) supports two types of low energy excitations. Excitations of the 
first type preserve the number of pairs (pair-preserving excitations). The second type of 
low lying excitations (pair-breaking excitations) is obtained by breaking a single electron 
pair. In the thermodynamical limit both types of excitations are gapped with the same gap, 
A p = A b = 2Ao, where A p and A b are the energy gaps for pair-preserving and pair-breaking 
excitations respectively. In Section we evaluate leading finite size corrections (of order 
1/m) to the gaps A p and A b . Interestingly, it turns out that these corrections coincide, even 
though the two gaps are not identical in higher orders in 1/m. In the limit A Q /D —>■ 0, our 



result yields A p = A b = 2Ao — d. We also show that the energy levels of lowest excitations 
of two types cross at certain value of the coupling constant A. 

Another measure of the low energy properties of BCS model (|l.lj) is the parity parameter- 
introduced by Matveev and Larkin. This parameter is defined as 



Aml = Eft* 1 ~ 7T (ElT 2 + El m s ) (1.5) 
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where E l s is the ground state energy of BCS Hamiltonian (jl.lj) with I electrons. Matveev 
and Larkin evaluated Aml in the physical limit Aq/D — > in two different regimes: Aq ^> d 
and Ao <C d. They found that in the first regime the leading finite size correction to the 
parity parameter (jl.5j) comes entirely from the stationary point (mean field) expression for 
the ground state energy of BCS Hamiltonian (|1.1JI . Here we use our method to calculate Aml 
in the regime A > d for an arbitrary ratio A Q /D. We show that the contribution of quantum 
fluctuations to the leading finite size correction to A M l behaves as (A /D) ln(A /D) for 
small Aq/D. 

The ground state energy of pairing Hamiltonian (jl.lj) has been discussed recently in a 
number of papers. Numerical fits for finite size corrections to the ground state energy in the 
weak coupling regime, A C 1, have been proposed^. Here we evaluate the leading finite 
size correction exactly and find a complete agreement with numerical results^ in the weak 
coupling regime. 

In Ref. [7|, authors studied the condensation energy, defined as the difference between the 
ground state energy and the expectation value of BCS Hamiltonian (jl.lj) in the Fermi ground 
state. This difference was calculated in second order perturbation theory in A and compared 
to BCS expression Eg^' s (X) —Eg^ s (0). The authors found that the two expressions become 
of the same order when d < A < V ' Dd and interpreted this as a new, "intermediate", 
regime of pairing correlations in metallic grains. We argue bellow that, although the finite 
size correction to the condensation energy indeed becomes of the same order as the BCS 
result for d < A < yDd, this fact does not indicate a new physical regime, but is rather 
an artifact of the model. Main contribution to the finite size correction to the condensation 
energy comes from energies close to the ultraviolet cutoff D and therefore is beyond limits 
of applicability of BCS Hamiltonian (jl.lj) . Effects coming from this range of energies can 
be properly accounted for— within the Eliashberg theory 14 . 

The paper is organized as follows. Section|H]is devoted to the review of a general method 10 



of 1/m expansion due to Richardson. In Section ITTT) we show that Richardson's results can 
be used to evaluate ground state and excitation energies of BCS Hamiltonian (jl.lj) to any 
order in 1/m and explicitly calculate the leading correction to the ground state energy. In 
Section I1VI we discuss various limits of our results and make a comparison with previous 
work. Results for the excitation spectrum and Matveev Larkin parameter are collected in 
Sections IS and IVII respectively, where we also determine the gaps for pair-breaking and 
pair-preserving excitations and discuss the range of applicability of the 1/m expansion. 

II. REVIEW OF RICHARDSON'S 1/m EXPANSION 

Here we briefly review Richardson's 1/m expansion 10 for the ground state and excitation 
energies of pairing Hamiltonian (jl.lj) . The details can be found in the original work 10 . 
In subsequent sections we will use Richardson's results to explicitly evaluate finite size 
corrections to the low energy spectrum of BCS Hamiltonian (jl.lj) . 

Richardson's 1/m expansion is based on an electrostatic analogy to equations (II. 2j) . In 
this analogy, the roots Ei of equations (jl.2j) are interpreted as locations of m two-dimensional 
free charges of unit strength in the complex plane. The free charges are subject to a uniform 
external field —l/(Xd) and the field of n fixed charges of strength 1/2 located at the points 
€k on the real axis. The total electrostatic field at a point z associated with the charge 
distribution is given by 

m 1 1 n 1 1 

^) = gi^- 5 £— - Yd (2.D 

The field F(z) contains complete information about the spectrum of BCS Hamiltonian 
(jl.lj) . For example, the energy spectrum is related to the quadrupole momentum of F(z). 
Indeed, defining multipole moments of F(z) by 

oo 

F(z) = J2 F {m) z- m (2.2) 

771=0 

and expanding equation (J2.1J) in 1/z, we obtain 

m n 

E = 2Y / E i = 2F^ + J2e k (2.3) 

j=i fc=i 

_ _L = F (o) (2 .4) 



m-- = F (1) (2.5) 

The 1/m expansion is facilitated by the following field equation that can be derived from 
equations (|1.2|) and (|2.1j) : 

dF | p^- 1 ^ ! | i/V ! | 2 V v Hk (26) 

dz^ 2^{z-e k )^ A\^z-e, Xd) ^ z - e k K ' } 

where H k is the field at the location of the fixed charge e& due to the free charges 

H k = E-\ (2-7) 

Equation (j2.fi j) can be solved by expanding the field F(z) in powers of 1/m 

oo 

F(z) = J2F r (z) (2.8) 

r=0 

where F r (z) is of order m 1 ~ r . It turns out 10 that the lowest order in 1)2.8)1 . Fq(z), together 
with field equation (J2.6J) completely determine the field F(z) to higher orders in 1/m. More- 
over, to obtain higher orders, F r (z) for r > 1, from .Po(z) one needs to solve only algebraic 
equations. 

Different states of the system are described by different F (z). For example, one can show 
that the BCS ground state corresponds to 



J(z - /i) 2 + A 2 

F (z) = - £ vl N (2.9) 

* 2(^-e fc ) v /( efc - / i) 2 + A 2 

The parameters A and /i correspond to the BCS gap and chemical potential respectively. 
Equations for A and \i can be derived by substituting F (z) into equations (|2.4|) and (|2.5|) 

^ = E / l (2- 10 ) 

n -2m = E ^"^ (2-11) 



V v /(e,-/i) 2 + A 2 

There are no higher order corrections to equations (J2.10J) and 1)2.11)1 . since by construction 
Fq(z) yields exact monopole and dipole moments of F(z), F^°'(z) and F^ l \z). 

Note that, according to equations 1)2.1)1 and (J2.9J1 . F (z) also describes the fixed charges 
exactly, since 

(2.12) 
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-e k )F {z) = 





Z~*e k 
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Higher order corrections to the field F(z) can be expressed only in terms of £&, A, /z and 
finite zeroes of F (z) 



E 



i 
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i(x,-e fe W(e fe -/i) 2 + A 2 



For example, 



Fx(z) 



1 /„ z + e k -2/j „ 2 + xi - 2[i z — \x 
2Z[z) [^ Z(z) + Z(e k ) " ^ Z(i) + Zfa) ~ ~zjz) \ 



(2.13) 



(2.14) 



where 



Z(z) 



[z - fi) 2 + A 2 

One can show (by e.g. sketching the LHS of equation (J2.13J) ) that there are n — 1 finite 
solutions to equation (J2.13J) . each of them lying between two consecutive single electron 
levels e k . 

The ground state energy to the first two orders in 1/m, i.e. to the order m°, can be 
obtained from Fq(z) and Fi(z) using equation ()2.3|) . 



E = E + E 1 

A 5 



^, = E e *-M"-2m) + — - J] ^ - /z) 2 + A 2 
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z=i 



(x,- Ai ) 2 + A 2 



Pi 



(2.15) 
(2.16) 



where 
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(xi - e k y 



PI=T 



r(x,- efe )V(e fc -/i) 2 + A 2 



To calculate excitation energies one needs to appropriately modify Fq(z), the lowest order 
in 1/m of the electrostatic field F(z). Here we simply write down excitation energies to the 
first two nonzero orders in 1/m referring readers interested in the detailed derivation to the 
original work- 
ed) = ei(/) + e 2 {l) l = l,...,n-l 

(2.17) 
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ei (/)=2^(x z -/i) 2 + A 2 



where 



F{(z) = F^z) 



(s, - /z) 2 + A 2 



(z-xO v /(^- ; u) 2 + A 2 z-zj 



(2.18) 



(2.19) 



and e(l) is the excitation energy relative to the ground state. 

Finally, we note that the lowest nonzero order of 1/m expansion, Eq and e\(l) for the 
ground state and excitation energies, reproduces the mean field (BCS) results for pairing 
Hamiltonian 1)1.1)1 . Therefore, the mean field for pairing Hamiltonian ()1.1|) is exact in the 
thermodynamical limit, while contributions E\ and e2<7), equations 1)2.16)1 and 1)2.18)1 . are 
leading finite size corrections to the thermodynamical limit. 

III. GROUND STATE ENERGY 

Here we evaluate the leading finite size correction to the ground state energy of BCS 
Hamiltonian 1)1.1)1 . 

First, we note that, as shown in Appendix A, expression (J2.16J1 for the finite size correction 
E\ can be cast into a simpler form 

/„ x n— 1 n 

Ex = Xd r - m + £ y/( Xl - ^ + A 2 - J2 y/(e k -»)* + A* (3.1) 

VZ J 1=1 k=X 

To facilitate comparison to the mean field BCS result (J1.4J1 . we assume bellow n = 2m equally 
spaced single electron levels e*, = (k — m — l/2)d with energies ranging from D = (m— l/2)d 
to —D. It should be emphasized, however, that explicit results in terms of A, /i, and the 
density of states u(e) can be equally well obtained for arbitrary continuous v{e). 

Since n = 2m and e^ are distributed symmetrically with respect to zero, equation (J2.11J) 
yields // = 0, while equations ()2.10)l . ()2.15)1 . and ()3.1)1 become 

9 2m -i 

Xd = ti y/el + A^ (3 ' 2) 

A s 



2m 



£ o = ^-EV^ + A2 (3-3) 
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«i=E^ + a 2 -E V4 + A2 (3-4) 

1=1 k=l 



Equation ()2.13)) for x\ now reads 



2m i 

/(*i) = E ; /2 A2 = (3.5) 



Since for each e^ there is ty = — e^, /(^) is an odd function of z. Therefore, x\ = is 
a solution of equation (J3.5J1 . while the remaining n — 2 = 2m — 2 nonzero solutions come 

8 



in pairs of x\ and —x\. Let us label m — 1 positive roots x\ with I = 1, . . . , m — 1 and 
relabel m positive single electron energies e/c with fc = 0, 1 . . . , m — 1. Then, we can rewrite 
equation (|3.4|) as 

m— 1 



£?i = A 



/rf 2 



m— 1 



(3.6) 



A2 + 2E V^ + A2 -E V e ' + A2 

Y ■* i=l fe=l 

where we have separated contributions to the summations of x/ = and e& = ±d/2. 

Because xi is located between q and q_i = ei — d, we can write it as xi = ei — aid, where 



< oci < 1. Expanding wa; 2 + A 2 in a^ in the vicinity of x\ = e\ and bearing in mind that 
d rj D/m is of order 1/m, we obtain 



m— 1 



£? x = -a - 2 x; 



(3.7) 



S Ve? + A 2 

where we neglected terms of order 1/m. With the same accuracy, we can replace the sum- 
mation over k with an integration 

cD 



E t = -A 



2/ de^L= 
Jo Ve 2 + A 2 



(3.8) 



Note that E\ is indeed of order m° as it should be. The function a(e) is evaluated in 
Appendix B. The result, up to terms of order 1/m, is 
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We 2 + A 2 - e^D 2 + A 2 



ZVe 2 + A 2 + e^D 2 + A 2 



(3.9) 



Introducing a new variable 



.r 



iln 

7T 



DVe 2 + A 2 - eV^D 2 + A 2 



.ZVe 2 + A 2 + eVD 2 + A 2 

we can cast expression (|3.8|) into a more convenient form 

i" 00 dx 



DA sinh(7rx/2) 
A 2 cosh 2 (vrx/2) + L) 2; 



(3.10) 
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(3.11) 



'A 2 + L> 2 (cosh(7rx/2))" 

To complete the evaluation of the ground state energy to order m°, we also need to calculate 
the leading term, E with the same accuracy. The first step is to replace summation in 
equations (|2.10|) and ()3.3|) with integrations according to the following formula: 

"■2 rri2d d 

f{jd) = / dxf(x) + 



J 12 r-n^d A 

d E f^d) = / dxf(x) + - [f(md) + fined)] + o(l/m) 



Equations (J2.1UJI and (j3.3|) now read 

A ~ J-D Ve 2 + A 2 + VA 2 + D 2 ^ ' ' 

E = ^---f D deVe 2 + A 2 - VA 2 + D 2 (3.13) 

Act a j-d 

The solution of equation (J3.12J) for A to order m° is obtained by dropping the second term 
on the RHS. Evaluating the integral, we obtain A = -D/[sinh(l/A)] in agreement with 
equation (|1.4j) . To compute the correction of order 1/m to A, we substitute A = Ao + SA 
into equation ()3.12|) and expand in SA. Keeping only terms of order 1/m, we find 

A = A + d^ (3.14) 



Plugging A into equation (J3.13)) and using J A 2 , + D 2 = D coth(l/A), we obtain up to terms 
of order 1/m 

E = -(m + -\Dcoth(l/\) (3.15) 

Note also that A in expression (|3.11|) for Ei can be replaced by A up to terms of order 
1/m. Thus, the ground state energy of BCS Hamiltonian (jl.lj) for m pairs and n = 1m 
equally spaced levels is 



E g ^ = -£>coth(l/A) 

where 



m + - + 0(A) 



(3.16) 



DC 



0(A) =2/ n dX ,^ - C ° sHnx/2) (3.17) 

J o ^l 1 + x ) ^/coshVx/2) + sinh 2 (l/A) 

The plot of function 0(A) is shown on Fig. 1. 

The finite size correction to the mean field BCS result (jl.4j) is 



E g . 8 . = E^ b + E f . s . E f . s . = -D coth(l/A) 



\+m 



(3.18) 



Note that Ef s is different from Ei given by equation (J3.11J1 due to contribution of order m° 
from i?o- 

Higher order corrections to the ground state energy can also be evaluated explicitly. The 
first step is to express them in terms of A and Xi following prescriptions of Ref. [lfj. Then, 
A and xi have to be calculated to appropriate order in 1/m using methods of this section 
and Appendix B. Final results for higher order corrections will involve multiple integrations 
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FIG. 1: The plot of function (ft(\) defined by equation l|3.17|) . This function appears in leading 
finite size corrections to ground state (|3.16|) and excitation (|5.1Ull5TT3*|) energies of BCS Hamiltonian 
(|1.1|) and to Matveev-Larkin parameter (|6.2j) . Note the asymptotics (j)(X) — > and 4>(X) — > 1 for 
A — » and A — » oo respectively. 

similar to the integration in equation (|3.17|l . For example, the expression for the correction 
of order 1/m contains a triple integral. 

The general case when the distribution of single electron levels in the limit m,n — > oo, 
m/n = fixed is described by a continuous density of states u(e) can be treated similarly. 
Final expressions for corrections will now be in terms of A, /i, and v(e). For example, the 
correction of order m° will be again given by the integral in equation (|3.8|) where the limits 
of integration should now be — D and D, e has to be replaced with e — /i, and the integrand 
has to be multiplied by u(e) . The function a(e) will still be given by equation (jB.2|) where 
now u(e) has to be included under the integral. 

IV. COMPARISON TO PREVIOUS STUDIES 

Here we analyze our result and compare it to previous results. First, we check whether 
equation (J3.18J) reproduces the results of 1/A expansion 15 around A = oo. Expanding 
the integrand in equation (|H.17jl in 1/A, evaluating the resulting integrals, and expanding 
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coth(l/A) in 1/A, we obtain 
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Comparing this expression with terms of order m° in 1/A expansion 15 for the ground state 
energy (see equation (30) of Ref. ll5T). we find that the two results coincide. 

Now let us consider the limit of small A. The asymptotic behavior of 0(A) for small A is 
worked out in Appendix C Here we write down the first two terms 

0(A) = A + ln2-A 2 + 0(A 3 ) (4.1) 



Expanding coth(l/A) = J\ + Aq/D 2 in Aq/D and using D = (m — l/2)d, we obtain from 
equation ()3.1fij) 

E g .,. = -D (m + -) - ^ - D\ - \n2 ■ D\ 2 + 0(A 3 ) (4.2) 

The first term in equation (J4.2J) is the energy of noninteracting Fermi ground state to order 
m°. The second term is the nonperturbative mean field (BCS) contribution to the ground 
state energy. The first two terms are extensive and survive the thermodynamical limit. The 
last two terms give the correction to the ground state energy that one would obtain in the 
second order of ordinary perturbation theory in A around noninteracting Fermi ground state. 

We see that our result ()3.18|) yields the leading finite size correction to the thermodynam- 
ical limit for all values of A. In particular, there is no breakdown in the regime of ultrasmall 
grains, i.e. for d > A . As we will see in subsequent sections, this is not a generic feature 
of our approach, but is specific to the ground state energy and is probably related to the 
ultraviolet nature (see bellow) of the finite size correction calculated above. 

A frequently discussed quantity 2,6 ' 7 is the difference between the ground state energy and 
the expectation value of BCS Hamiltonian (jl.ljl in the unperturbed Fermi ground state, 
\F.g.s.), i.e. a state where single particle levels bellow the Fermi level, e& < 0, are doubly 
occupied, while the remaining levels are empty. This difference is often called condensation 
energy, even though this name is misleading for the reasons detailed bellow. However, to 
facilitate a comparison with results of Ref. |y and [J we will use the same terminology in this 
section. We have 

E cond . = (F.g.s.\H BC s\F.g.s.) - E g . s . = -D [m + - J - 2Am<i - E g . s . 

12 



Using D = (m — l/2)d and equation (|4.2jh we obtain 

^cond. = ^ + In 2 • DA 2 + 0(A 3 ) (4.3) 

Comparison shows that the exact result (J4.3J) for E con d. to order m° is in complete agreement 
with fits to numerical data££ 

Finally, note that the second term in expression (J4.3J) is ultraviolet divergent, since it 
depends explicitly on the ultraviolet cutoff D. For pairing by phonons the ultraviolet cutoff 
D can be identified with the Debye energy ujd- To properly take into account any effect 
that comes from energies of the order of up, one needs to go beyond the BCS theory which 
is appropriate only at energies much lower than ujd- The contribution from these energies 
to finite size corrections can be adequately treated 13 within the Eliashberg theory^- In 
particular, the hard cutoff at D = u>d has to be replaced by a soft effective cutoff due to the 
1/uj 2 decay of the phonon propagator for frequencies to ^> ujd- Therefore, even though the 
contribution of the finite size correction in equation (J4.3J) becomes important for Ao < yDd, 
the conclusion of Ref. |7| that this is an indication of any new physical regime is not justified. 

V. EXCITATION ENERGIES 

In this section we evaluate leading finite size corrections to lowest excitation energies. As 
we will see bellow, the results of this section are accurate only in the regime of relatively 
large grains, A > d, i.e. within terms of order o(d/A ). These higher order corrections can 
also be straightforwardly calculated using methods of Sections IIHI However, we will only 
evaluate corrections of order d/A here. 

As in Section IIHI we will perform calculations for the case of 2m electrons and n = 2m 
equally spaced levels €& = (k — m — l/2)d with energies ranging from D = (m — l/2)d to 
— D. In this case, equation ()2.11|) implies \x = 0. A more general case when the single 
electron levels are distributed with a smooth density of states can be treated similarly (see 
the discussion bellow equation (J3.I8J) ). 

Note that Hamiltonian (jl.lj) conserves the number of paired electrons. Therefore, the 
excitations can be grouped into two types: those that preserve the number of pairs and those 
that break pairs. Energies of low lying pair-preserving excitations in the thermodynamical 



13 



limit are given by equation (|2.17|) with \i = 



e[ = 2^xf + A 2 (5.1) 

where Xi are the roots of equation (|2.13|) . Low lying pair-breaking excitations are obtained 
by breaking a single pair and placing the two unpaired electrons on two single electron levels 
e a and e&. The energy of this excitation according to equation (J1.3J) is 

e b = e a + t b + E g . s . (e a , e b ) - E g . s . (5.2) 

where E g . s (e ai e b ) is the ground state energy of BCS Hamiltonian (jl.lj) with levels e a and e b 
blocked. In the thermodynamical limit, using equation (J2.15J) . we obtain 



el + A 2 + sjel + A 2 (5.3) 

Therefore, in the thermodynamical limit both types of excitations are gapped with the same 
gap 2A , i.e. 

A? = A? = 2A (5.4) 

Since pair-breaking excitations are capable of carrying spin-1, A 6 can also be called the spin 
gap. To calculate corrections to Af and A\, one needs to go beyond mean field approxima- 
tion. 

First, let us determine the energy of lowest lying pair-breaking excitations to order 1/m. 
Breaking a pair changes both the number of pairs to m! — m — 1 and also the number of 
unblocked levels to n' = 2m — 2 = 2m! . The lowest energy is archived by blocking levels 
£a — d/2 and e b = —d/2. Since this leaves the distribution of single particle levels symmetric 
with respect to zero, the chemical potential /i in equation (J2.11J) remains equal to zero, 
// = /i = 0. However, the blocking affects the gap A', since now terms corresponding to 
€k = ±d/2 have to be excluded from gap equation (J2.10J) . Using equation (j2.1()|) . we obtain 

1 2 1 

E / = / +E / ( 5 - 5 ) 

* fl + A' 2 VW4 + A2 V f\ + A^ 

where A' is the value of the gap with levels ±d/2 blocked. Expanding the LHS of equation 
(15. 5|) in 5 A = A' — A and using gap equation (j2.1()|) . we obtain 



6A = -dJl + ^ (5.6) 
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According to equation (J5.2J1 . to order 1/m the lowest lying pair-breaking excitations have 
the following energy: 

A b = E' {A') - E (A) + ^(A',xi) - £i(A,s,) (5.7) 

where -Eo(A) and i?i(A, xj) are given by equations (|3.3jl and (|3.4|) respectively and primes 
denote quantities for the ground state with levels ±d/2 blocked. Equations (|3.3JI . Q3.4J) . 
(loTfl) . and (l3~TiJ) imply 



£ ;(A') - B„(A) = 2A' + £ 4(e / + A f 2)3/2 = 2A + ^ - c^l + § (5.8) 

where -Ei(A) is given by equation (|3.11j) and <5x; is the change in x\ due to blocking levels 
±d/2. 

We see from equation (|3.5|) that the effect of removing levels e k = ±d/2 from the sum- 
mation in equation (J2.13J) is strongest for the roots closest to the blocked levels ±d/2. For 
these roots Sxi ~ d. On the other hand, due to an additional factor of x\ in front of Sxi 
in equation (|5.9|h the contribution of each of these x\ to the RHS of equation (J5.9J) is of 
order d 2 / A. By splitting the sum in equation (J3.5J) into two sums as in Appendix B, one 
can show that the contribution of all these roots to the sum in equation (|5.9JI is of order 
o{l/m). For the remaining roots, Sxi/xi is of order 1/m and each term in equation (J3.5J1 
can be expanded into Sxi/(xi — e&). We have 

i 2m i o 2m -| 



-^— ' r~o '. r~^/ , s ^-^ r~o : rrr/ <. s ™_ a ^-^ 



e ^±d/2 ^ + A 2 (xJ - e fc ) *=i y/? k + A\x h + <5x z - e k ) ^ A fc=1 ^2 + A 2( X/ _ ejfe) 
Expanding into Sxi, we obtain 

1 2 



<^E 



/ef, + A 2 (x, - e k y xA 

The summation here can be evaluated in the same way as the first sum in equation (JB.1JI 
of Appendix B is evaluated. Recall that roots of equation (J3.5J) Xi and therefore 5xi are 
distributed symmetrically with respect to zero. Using the notation introduced in the text 
following equations (|3.5|) and (J3.6J) . we have for xi > 



2d 2 J^ + A 2 sin 2 vr«(eO 
OX1X1 = ——. 

A 7T 2 
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where a(e{) is given by equation (|3.9|) . Substituting 5xiXi into equation ()5.9|) and using 
equations ()5.8|) . (J5.7J) . (J3.9J) . and ()3.11|) . we obtain 



/ A 2 <7A 
A 6 = 2A - d^l + ^ + -^ [1 + 0(A)] (5.10) 

where we used the change of variables (J3.10J) and 0(A) is defined by equation ()3.17|) . Expres- 
sion ()5.10j) yields the energy of lowest lying pair-breaking excitations up to terms of order 
o(d/(min[D,A ])). 

In the physical limit of weak coupling, A /D —* 0, according to equation (|4.1jl . expression 
(I5.10|) becomes 

A 6 = 2A - d + o(d/A ) (5.11) 

Next, we turn to excitations that preserve the number of pairs. Energies of these ex- 
citations to order 1/m are given by equations (|2.17|) and (|2.18|) . Equation (J5.1J1 shows 
that the lowest lying excitation corresponds to xi = 0. We have, up to terms of order 
o(d/(min[D,A ])) 



A p = 2A + 2 Y, 4 



P/ 



(F x ') 2 - (F) 2 + -^ - iV ' 2F{ 



dz x. 



in 



(5.12) 



where Fi(z) and F x '(z) are defined by equations (J2.14J) and (J2.19J) . Taking into account that 
both €k and x\ are distributed symmetrically with respect to zero and fi = 0, we can rewrite 
these equations as 



2Vz 2 + A 2 \y y/z 2 + A 2 + J e 2 k + A 2 V v^ 2 + A 2 + Jx 2 + A 2 V* 2 + A 2 



Fi(z) = F 1 (z) 



xf + A 2 l 



(z-a:i)v / ?TA I 2-xz 
Summations in Fi(z) and in equation (|5.12|) can be evaluated in the same way as sums in 
equations (J5.9J) and (J3.6J) have been evaluated. Even though this calculation looks rather 
different from the one that lead to equation (|5.1()jl . it yields an identical result, i.e. 

A p = A fe + o(d/(min[Z),Ao])) (5-13) 

Thus, both gaps coincide up to terms of order o(l/m). However, this coincidence is not 
preserved in higher orders. Indeed, it was shown in Ref. [l5j that in the strong coupling 
limit, A ^> 1, the gap for pair-breaking excitations is larger A 6 — A p ~ d 2 /A > 0. On 
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the other hand, at A = the gap for pair-preserving excitations is larger, A b — A p = —d. 
Therefore, the lowest energy levels of the two types of excitations cross at certain value of 
A . Equation (|5.13|) shows that the distance between the two levels is reduced from d at A 
to o(d/A ) d even when d -C A -C D. However, the knowledge of higher order corrections 
to the gaps A b and A p is needed to determine whether the crossing occurs in the physical 
regime A /D — > 0, i.e. at A — d. 

VI. MATVEEV-LARKIN PARAMETER 

Finally, let us evaluate the Matveev-Larkin parameter—. This parameter is a measure of 
a parity effect in the grain and is defined as follows: 

A ML = ^- +1 - \ {E 2 9 T 2 + Bfc) (6.1) 

where E l s is the ground state energy of BCS Hamiltonian ()1.1|) with I electrons. 

The calculation of Aml is similar to the one that lead to equation (|5.10|h only now we 
also have to take into account the change in the chemical potential 



/ A^ 
H>2m+2 — /J>2m = 2(/i2m+l — /^2m) = — 2(A2 m +2 — A2 m ) = d\J 1 + jS 



A 2m+2 - A 2m = O(d 2 /A ) 
The calculation results in 



A 6 . d A^ dA ,.,.. 

A ML = _ = A --^l + -^ + ^[l + 0(A)] (6.2) 

where 0(A) is defined by equation (|3.17|) . As before, this expression is accurate up to terms 
of order o(d/(mm[Ao,D])). In the physical limit AnjD — > 0, expression (|fi.2|) . according to 
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equation (|4.1jl . reduces to the one obtained in Ref. 

Aml = A - - + o(d/A ) (6.3) 

The first three terms on the RHS of equation (|6.2jl come from the mean field (stationary 
point) approximation (J2.15J) for the ground state energy. The last term in equation ()6.2|) 
represents the contribution of order 1/m of quantum fluctuations around the stationary 
point. The asymptotic behavior of this term in the physical limit A /D — > is given by 
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equation (|4.1|) . In terms of d, Aq, and D it reads d\n(Ao/D)Ao/D. In this limit quantum 
fluctuations will contribute to higher orders in d/ Aq as evidenced by the result^ for Aml 
in the regime d -C A . Therefore, it is of certain interest to use methods of Section ITTT1 to 
evaluate further corrections to A ML . 

We conclude this section with a comment on the range of applicability of 1/ra expansion 
detailed in this paper. It is clear from equations (J5.11J) and (J6.3I) that the expansion is 



applicable in the regime Aq > d. In fact, results of Ref. 



15l and [l0| (see also Section ITT]) 



suggest that the expansion is in powers of d/A with a convergence radius d/A Q ~ 1. 

VII. CONCLUSION 

In this paper we have shown that finite size corrections to the thermodynamical limit for 
pairing Hamiltonian (jl.lj) can be evaluated explicitly in terms of the BCS gap A , chemical 
potential /i, mean level spacing d, ultraviolet cutoff D, and the thermodynamic density of 
states u(e) to any order in d/A ~ 1/ra. 

We evaluated leading corrections to the ground state and lowest excitation energies, and 
to Matveev-Larkin parameter (equations (J3.181 15.101 15.111 15.131 16.21 I6.3J1 ). Our results for 
the ground state energy are in agreement with previous numerical studies. We showed that 
the finite size correction to the condensation energy is ultraviolet divergent and therefore 
comparing it to the BCS result is not justified. 

We found that the gaps for pair-breaking and pair-conserving excitations of pairing Hamil- 
tonian (jl.lj) coincide up to terms of order o(l/ra), where m is the number of electron pairs 
on the grain. In higher orders in 1/m the two gaps are different, the difference being of 
order d 2 /A , where d is the mean level spacing and A is the BCS gap (J1.4JI . We showed 
that the energy levels of the lowest excitations of two types cross at a certain value of the 
coupling constant A. 

The range of applicability of 1/ra expansion detailed in the present paper is A > d. In 
fact, we believe that in the physical limit A /D — > the expansion is a power series in d/A 
with a convergence radius of order one. 

Note that our results significantly simplify in the physical limit Aq/D —> (e.g. compare 
equations (J5.1UJI and (J5.11J) ). An interesting open problem is to take this limit directly in 
Richardson's equations (jl.2j) and to develop a simplified version of the 1/m expansion for 
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this case. In particular, this might help to address the problem of the crossover between the 
fluctuation dominated (d ^> A ) and the bulk (d <C A ) regimes. 
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IX. APPENDIX A 

Here we show that expression (J2.16|) for the correction to the ground state energy can be 
simplified to equation (|3.1|) . Indeed, define 

n d 1 

f( z ) = X! where d k 



Equation (j2.13|) now reads f{x{) = 0. The function f(z) has n — 1 finite zeroes at z = x\ 
and also a zero at z = oo. Its dual function, g(z) = l/f(z), has n — 1 poles at z = x\ and 
also a pole at z = oo with a residue (J2k=i ^fc) _1 - Therefore, it can be represented as 

( ) - v mi - - V mi +■ — 

f^z-xi £ fe 4 f^z-xi 2 

where we have used J2kdk = 2/ (Ad) in accordance with gap equation (J2.10J1 . The following 
equations for the residues of g(z) and f(z) are helpful: 



mi = lim \{z — xi)g(z)} = lim 



*-*i L *-»*« f( z ) f'(xi) 






1 



P 



1 1 mi , Ad 



4 lirn^J^ - e fe )/(z)] , (xj - e k ) 2 2 

where the prime denotes the derivative with respect to z. Using these equations, we obtain 



n— 1 at n—1 n 



Finally, substituting equation (|A.lj) into expression (J2.16J1 . we obtain equation (|3.1|) . 



19 



X. APPENDIX B 



In this Appendix we solve equation (J3.5J) for xi. As was discussed bellow equation (|2.14|) . 
each solution x\ lies between two consecutive single electron levels e&. Consider the solution 
x(e) that lies between e — d and e, where we dropped subscripts for simplicity. 

Now let us multiply equation (|3.5|) by d and rewrite it as 

d ^^ d 



E 



+ E 



o 



(B.l) 



\e k -e\<jd (x(e) - e fe )\/ e fc + A2 k h -e|>Jd (»(e) - ^k)\J4 + A 2 
where 1 € J € A = min[A,D]/<i. For example, one can choose J = yA. In the first 



summation in equation (jB.l|) . we 2 , + A 2 can be replaced by \/e 2 + A 2 with a relative error 
of order Jd/A. We obtain 



E 



d 



Kfc- 



-e|<Jd (xle 



e fc Wef + A 2 



-°x 



1 J 



p+l-a(e) p + a(e) 



where a(e) is defined by x(e) = e — a(e)<i. To determine a(e) to the leading (m°) order 
in 1/m, we can now take the limit m — ► oo. With a suitable choice of J (e.g. J = vA), 
J — > oo and (Jd)/D —> in this limit, while the second sum in equation (|B.1|) becomes a 
principal value integral. Using, 

1 1 



E 

p=0 



p+1 — a(e) p + a(e) 



we obtain 



7rcot(7ra(e 



r D 



— 7r cot(7ra(e)) 



de' 



(B.2) 



'-D (e-e'Ve^ + A 2 
Finally, evaluating the integral, we arrive at equation (|H.9jl . 

Corrections <5a(e) to a(e) of order 1/m and higher can also be evaluated explicitly by 
expanding equation ()3.5|) in 5a(e). These corrections contribute to terms of order 1/m and 
higher in the ground state energy. 



XI. APPENDIX C 



Here we determine the asymptotic behavior for small A of the integral 

dx cosh(7nc/2) 



0(A) 



< l + x2 ) ^cosh 2 (7nr/2) + sinh 2 (l/A) 
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First, we note that up to terms of order e 1 / x , one can rewrite this integral as 



oo 



dx 2 

0(A) = 2 / — — = where x = — r 

' 7T[l + (x + Xo) 2 ]v / TTF^ VTA 



oo 



Let us divide the domain of integration into three intervals: (—00, —a), (—a, a), and (a, 00), 
where 1 < a C io, and denote the corresponding integrals by I 3 , I 2 , and I\ respectively. 
Each of the integrals Ik can be expanded into its own small parameter that depends on a. 
The dependence on a will cancel out when the results are added together. We have 

/•oo r]rn poo rl-r p~ 7TX /'2 P ~ 3irx/2 _i_ 

1=2 = 2 — - = 0{e^ 7Ta/2 ) 

3 Ja 7T[1 + (x - X ) 2 ]Vl + e™ A 7T 1 + (x - X ) 2 ^ 

* c?x 2 /- a dx 

a 7r[l + (x + x) 2 ]^! + e-* x x 2 J -a 7Ta/1 + e~ nx 

xdx 27ra + 41n2 ^ / a 2 \ 
+ ... = ^ + 0' 



XJJ 7-a 7T\/l + C~ nX n2x \ x 



:x 



dx r™ dx 1 - e~™/2 + .. . 2 2a /a 



r __ o / = o / / _i_ /n 

ia 7r[l + (x + X ) 2 ]\/l + e"' 1 ia 7T 1 + (x + X ) 2 7TX 7TXq V^O^ 

Adding Ji, J 2 , and I3, we obtain equation ()4.2|) . Higher order terms can also be calculated 
by the same method. 
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